Synchronization in a neuronal feedback loop 
through asymmetric temporal delays 
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We consider the effect of asymmetric temporal delays in a system of two coupled Hopfield neurons. 
For couplings of opposite signs, a limit cycle emerges via a supercritical Hopf bifurcation when the 
sum of the delays reaches a critical value. We show that the angular frequency of the limit cycle is 
independent of an asymmetry in the delays. However, the delay asymmetry determines the phase 
difference between the periodic activities of the two components. Specifically, when the connection 
with negative coupling has a delay much larger than the delay for the positive coupling, the system 
approaches in-phase synchrony between the two components. Employing variational perturbation 
theory (VPT), we achieve an approximate analytical evaluation of the phase shift, in good agreement 
with numerical results. 
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Synchronization phenomena among coupled systems 
are abundant in nature [l|, [1J . The coupling is often not 
instantaneous; rather finite time delays exist. In gen- 
eral, time delays can cause an otherwise stable system 
to oscillate [!, 0, Q and may lead to bifurcation scenar- 
ios resulting in chaotic dynamics [1, 0]. For example, 
delay-induced oscillations have been reported for neural 
networks genetic regulatory networks @, and mod- 
els of population dynamics (Io| to name just a few. The 
delays for the different coupling mechanisms in such net- 
works do not need to be uniform, which may have an 
important effect on the system dynamics. For instance, 
it has been shown that distributed delays can stabilize 
a dynamical system (Tlj . and the influence of delayed 
inhibitory feedback has been studied In regard to 

network synchrony, the question arises under what condi- 
tions this special form of network behavior can be main- 
tained when the temporal delays are nonuniform. 

Asymmetric time delays in the visual pathway can be a 
pathological condition, as they are associated with many 
diseases However, when feedback loops in biologi- 

cal systems have evolved to feature different latencies for 
feed-forward and feedback projections, this might pro- 
vide a hint that asymmetric delays can also be beneficial 
to a system's functioning. In the avian visual system, 
the optic tectum is reciprocally coupled with the nucleus 
pars parvocellularis (Ipc), a subnucleus of the nucleus 
isthmi 14]. The coupled systems, tectum and Ipc, re- 
spond with synchronized oscillatory bursts to visual stim- 
ulation (l~5T |. Remarkably, the Ipc axons projecting to 
the tectum are thick and myelinated (fast action poten- 
tial propagation), whereas tectal neurons projecting to 
the Ipc possess comparatively thin axons and are un- 
myelinated (slow action potential propagation) [14j . The 
Ipc-to-tectum delay may thus be as short as a fraction 
of a millisecond, whereas the delay for the tectum-to-Ipc 
projection can be expected to be of the order of tens 
of milliseconds. It therefore seems natural to conjecture 
that the asymmetry in the delays may play a functional 



role in the feedback system. 

To explore this conjecture we investigate a model sys- 
tem of two coupled Hopfield neurons [161 ] with asymmet- 
ric delays, described by the coupled first-order delay dif- 
ferential equations (DDE's) 



du\(t) 

dt 
du 2 (t) 
dt 



= -Ui (t) + Cii tanh[u 2 (t — r 2 )] , 

= -U2(t) + a-z tanh[ui(t - n)} . (1) 



Here, u\ and m 2 denote the voltages of the Hopfield neu- 
rons and Ti and t 2 are the signal propagation or pro- 
cessing time delays, while a\ and a 2 describe the cou- 
plings between the two neurons. The system of DDE's 
(UJ) has a trivial stationary point at u\ = it 2 = 0, the 
stability of which has been analyzed in detail, e.g., in 
ref. [17| • For aia 2 < —1 the fixed point at the ori- 
gin is asymptotically stable as long as the mean of the 
time delays r = (n + 72)/2 does not exceed the crit- 
ical value To = sin -1 [— 2wo/(aia 2 )]/(2wo), where uiq = 
i/|aia 2 | — 1. When the sum of the delays is increased, 
the origin becomes unstable and a limit cycle emerges via 
a supercritical Hopf bifurcation at r = to. Note that the 
characteristic equation for the system (fT]), which deter- 
mines the condition for a periodic solution to exist, only 
depends on the sum of the two delays. A linear stability 
analysis can thus provide no insight toward a possible 
role of asymmetry in the delays. Furthermore, standard 
methods for bifurcation analysis, as described, e.g., in 
refs. [IjdH are only suitable for examining the nonlinear 
dynamical system in the immediate neighborhood of the 
bifurcation. In contrast to that, in this letter we aim at 
obtaining results that also hold for large delays, i.e., far 
away from the bifurcation. 

We first investigate the effect of asymetric time delays 
through numerical simulations. For a subsequent per- 
turbation expansion we define the expansion parameter 
e = V T ~ T o- Figure [T] shows numerical solutions of the 
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FIG. 1: Numerical solutions to the system of DDE's |T} for 
the choice of parameters a\ = — 1 and a 2 = 2 and for different 
values of the time delays ti, T2 (transients not shown). So- 
lutions for the case n = T2 are represented by solid lines. 
Dashed and dotted lines represent solutions for the cases 
Ti = and T2 = 0, respectively. For each set of lines the 
value of the delay parameter e = \/t — to increases from the 
innermost limit cycle (e = 0.1) to the outermost limit cycle 
(e = 1.0) in increments of Ae = 0.1. 
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FIG. 2: (Color) Plot of the phase shift between ui(t) and 
U2 (t). Numerical results for the scalar product <f> as given 
by ((2)1 are color coded for combinations of n and n with 
< T1/2 < 2.8. Red and blue indicate negative and positive 



1/2 

values of <j>, respectively, 
solution exists. 



In the black region, no periodic 



system of DDE's (JXJ) for different values of the time de- 
lays Ti and T2 and for the choice of parameters a\ = — 1, 
a 2 = 2. The amplitude of the limit cycle is only deter- 
mined by the value of e and thus remains unchanged when 
the temporal delays are chosen to be different. However, 
we observe that the phase between the periodic activi- 
ties of ui(t) and 112(f) does depend on the asymmetry of 
the delays. In order to quantify this phase difference, we 
consider the normalized scalar product 



j-Tq+T ^ U\{t)u 2 {t) 



fi° +T dt ui(t)ui(i) fr° +T dt u 2 {t)u 2 {t) 



1/2 



(2) 



Numerical results for this quantity are shown in fig. [2] 
We find that for time delays which are equal or at least 
not too asymmetric the scalar product <fi is approximately 
zero, which corresponds to a phase shift of 7r/2 between 
iti(f) and u 2 (t), assuming that they can be described by 
sinusoidal functions. However, when the delays are asym- 
metric, the scalar product <f> becomes larger in magni- 
tude, being negative for t\ > t 2 and positive for t 2 > t\. 
Specifically, for ti = the scalar product approaches 
unity for a growing delay t 2 , corresponding to in-phase 
synchronization between the the two components. 

We now aim at achieving an approximate analytical 
calculation of <fi. To this end, we first derive the per- 
turbation scries for the periodic solution u(f) and its 
angular frequency lo of the system ([1]) by applying the 
Poincare-Lindstedt method [l9j. Since a supercritical 
Hopf bifurcation occurs at r = To, we assume that 
the amplitude and frequency of the new periodic states 



are analytic in e and expand them as u(f) = eU(f) = 
e [U<°) (f) + eU«(f) + ...], w(e) = u; + ew x + e 2 co 2 + . . .. 
Furthermore, for convenience we introduce the rescaled 
independent variable £ = w(e)f and write U(f) = V(£). 
The expansion thenproceeds in a way very similar to the 
approach in ref. [201 ]. where the frequency of the limit 
cycle is calculated perturbatively for increasing mean of 
time delays. However, we introduce an additional param- 
eter fx, which is defined as the Ti-value of the intersection 
point in the t\-t 2 plane between the line that marks the 
boundary between the regions in which a periodic solu- 
tion does or does not exist, and a line perpendicular to 
this boundary through a given point (n, T2). To nth 
order in e, we have to solve a system of differential equa- 
tions of the form 



di 



dV 2 (n) (0 



to 



(3) 



+ ^y 2 ( " ) K~co(2T -n)] + / 1 (n) (e), 

UJQ UJQ 



where the inhomogeneity f ^ (£) is determined by the 
solutions to previous orders. Since we require that the 
solution V( n ) (£) be periodic in £ with period 27T, we can 
impose certain conditions on the inhomogeneity f 
Namely, we demand that f not contain terms that 
would lead to non-periodic solutions for V(™)(£), i.e., 
f must not contain secular terms. These condi- 
tions, which can be derived by expanding both the nth 
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FIG. 3: (Color) Perturbative results for the phase shift between ui(i) and U2(t). The color-coded plots show the difference 
between the numerical result from fig. [2] and the perturbative results up to order g 3 . 



order limit cycle solution V'™)(£) and the inhomogeneity 
f (")(£) into a Fourier series, read 

a 2 sin(o;oTo)ai™ 1 ) + sin[w (To - h)] 

+ (3^1 cos[w (t - fi)] = , 



» 



l { cos[cj (t - n)] - (3^1 sin[cj (T - h)] 

- a 2 sin^oro)^ = 0. (4) 

Here aj^ t and i denote the coefficients of the cosine 
and sine terms in the Fourier expansion of the inhomo- 
geneity f i respectively. Imposing these conditions 
on the inhomogeneity in ([3]) allows us to determine the 
angular frequency correction uj n and the Fourier expan- 
sion coefficients for V(™~ 2 )(£). To second order in e wc 
find 



ui 2 = 



ujqTq + cos(w r ) sin(w To) 



(5) 



while u>i vanishes. This value is identical to the one found 
in ref. [20( depending only on uiq and To but not on T\ 
or T2- Since this observation holds to all orders, we thus 
find that the period of the oscillations is independent of 
any asymmetry in the time delays. Furthermore, we find 
that only even perturbative orders lead to nonvanishing 
contributions for both the angular frequency w and the 
limit cycle V(£); we therefore define the new expansion 
parameter g = e 2 . Denoting the expansion to order g 



N 



of the quantity © by </> (Ar) , we find 



kW - 



cos[^o(2t - fi)] - ^osin[cj (2ro - fj)] 
sign(ai) v /l + u 2 



(6) 



-g- 



2(t -ti)uq {sin[w (2r - n )] + w cos[w (2r - f\ )] } 



sig 



;n(ai)[sin(2w To)+ 2a; r ](l + uj 2 ) 3 / 2 



Focussing on the choice of parameters a\ = — 1, a 2 = 2, 
which leads to wo — 1, T o — 7r /4, we can determine 
the expansion coefficients for (j)^ up to the third or- 
der. Figure [3] shows a comparison of our perturbative 



results and the numerical result from fig. [21 For small 
time delays, the accuracy of the results from the pertur- 
bation expansion is good and improves with increasing 
order. However, as g increases, the perturbative results 
cease to converge and no longer provide an acceptable 
approximation. As is typical for perturbative methods, 
our approach has yielded a divergent series. In order 
to improve the quality of our results, we now perform 
a resummation of the perturbative expansion employing 
variational perturbation theory (VPT). 

VPT is a nonperturbative approximation scheme based 
on a variational approach due to Feynman and Klein- 
ert [2l[, which has been systematically developed over 
the last few years, establishing its applicability in var- 
ious fields of physics [IS SI SI 0, SE Si- VPT 
permits the evaluation of a divergent series of the form 
f^Hg) — J2n=o an 9 n ana - yields a strong-coupling ex- 



m=0 "mi 



pansion of the generic form f(g) = g p / q ^2 
Here, p and q are real growth parameters characterizing 
the strong-coupling behavior. The convergence of the se- 
ries after resummation is exponentially fast and uniform 
with respect to other system parameters such as temper- 
ature, coupling constants, spatial dimensions, etc. [27| . 

In order to perform the resummation, one introduces 
a variational parameter K for the perturbation series ac- 
cording to Kleinert's square- root trick [23[. The series is 
thus transformed to the expression 



(7) 



fe=0 



k J \ K 2 



derived in detail in ref. [20]. The influence of the vari- 
ational parameter K is then optimized according to the 
principle of minimal sensitivity [28[; i.e., the optimized 
value is determined by solving for the roots of the 

first or higher derivatives of f^ N \g, K) with respect to 
K . The iVth order VPT approximation is then obtained 
by evaluating {3 at this optimized value: fyp^(g) — 
f( N ^(g, K( N >). This variational result generally holds for 
all values of the coupling constant g. Furthermore, by 
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FIG. 4: (Color) VPT results for the phase shift between u\(t") and uz(t). The first color-coded plot shows the first-order- VPT 
result as given by 1(5). The three other plots show the difference between the numerical result and the results from the first 
three orders in VPT. For clarity, this difference has been augmented by a factor of 10 and 100 in the results for the first and 
for both the second and third order, respectively 



considering the limit of large g, it allows the extraction 
of the strong-coupling coefficients b m . 

In our case of the perturbation series for </>, the values 
of the growth parameters p and q turn out to be the same 
as those that we determined in ref. [2(| for the angular 
frequency, namely p = — 2, q = 2. Our first-order result 
after resummation then reads 

(1) (2 + 7r)(l-2cosfisinfi) 
<Pvpt(3) = ( 8 ) 

x [(2+7r)(cosfi — sinfi) + g(ir — 4fi)](cosfi + sinfi)] -1 . 

The first color-coded plot in fig. [4] shows a graphical rep- 
resentation of this result. The agreement with the numer- 
ical result from fig. [5] is excellent. While the second VPT 
order provides a significant improvement when compared 
with the first order result, third order results are slightly 
superior to those of second order. 

In conclusion, our investigation of a neuronal model 
system shows that asymmetric temporal delays can con- 



trol the phase in a feedback loop and lead to synchronous 
oscillations. Specifically, in-phase and anti-phase syn- 
chrony arises when the delays are maximally asymmet- 
ric. Furthermore, after a variational resummation of the 
perturbation series for 0, we have a very accurate ap- 
proximate result for this quantity even in low orders and 
throughout the full t\-ti plane. 
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